#### Run the list experiment data
#### Assumes you're in the base folder of the replication file

library(tidyverse)
library(magrittr)

#dotwhisker has been archived, you may need to download an archived version of it.

library(dotwhisker)
library(broom)
library(stargazer)

### Grab data
dat <- read.csv("data/Election in Nepal  v1_3.csv")
polpart.fit <- readRDS("analysis/FApolpart-fit.rds")
wealth.fit <- readRDS("analysis/IRTwealth-fit.rds")

### List experiment
le <- dat[,c("group.1.clientism.measures.9.cli.mea.2.le.in", 
             "group.1.clientism.measures.9.cli.mea.2.cli.mea.list.experiment.count",
             "enumerator.entry.entry.region",
             "group.1.demo.5.dem.var.gender",
             "group.1.demo.5.dem.var.age",
             "group.1.demo.5.caste.sub.dem.var.caste",
             "group.1.demo.5.dem.var.religion",
             "group.1.demo.5.dem.var.readletter",
             "group.1.demo.5.dem.var.writeletter",
             "group.1.demo.5.dem.var.education",
             "group.1.socio.economic.7.soc.sta.house.material",
             "group.1.socio.economic.7.soc.sta.residence",
             "group.1.socio.economic.7.soc.sta.cook.source",
             "experiment.debriefing.10.exp.deb.qn.exp.deb.wealth",
             "experiment.debriefing.10.exp.deb.qn.exp.deb.guess.research.team")]
names(le) <- c("le.group", "le.response", "region", "gender", "age", "caste",
               "religion", "readletter", "writeletter", "edu",
               "house.material", "rural.urban", "cook", "self.relwealth",
               "who.conduct")

le$le.group <- as.factor(le$le.group)
le$any <- le$le.group != "1"

le$polpart <- apply(polpart.fit[[2]], 2, mean) ## posterior average
le$polpart.high <- le$polpart > mean(le$polpart)

le$wealth <- apply(wealth.fit$theta, 2, mean) ## posterior average
le$wealth.high <- le$wealth > mean(le$wealth)

#rough and quick NA replace for character variables since na_if package change
le$region %<>% na_if("dk")
le$region %<>% na_if("dta")
le$gender %<>% na_if("dk")
le$gender %<>% na_if("dta")
le$caste %<>% na_if("dk")
le$caste %<>% na_if("dta")
le$religion %<>% na_if("dta")
le$religion %<>% na_if("dk")
le$readletter %<>% na_if("dk")
le$readletter %<>% na_if("dta")
le$writeletter %<>% na_if("dk")
le$writeletter %<>% na_if("dta")
le$edu %<>% na_if("dk")
le$edu %<>% na_if("dta")
le$house.material %<>% na_if("dk")
le$house.material %<>% na_if("dta")
le$rural.urban %<>% na_if("dk")
le$rural.urban %<>% na_if("dta")
le$cook %<>% na_if("dk")
le$cook %<>% na_if("dta")
le$self.relwealth %<>% na_if("dk")
le$self.relwealth %<>% na_if("dta")
le$who.conduct %<>% na_if("dk")
le$who.conduct %<>% na_if("dta")

### some recoding
le %<>% mutate(edu2=(edu=="secondary-education" | edu=="above-intermediate-education" | edu=="intermediate-education"))

## Truncation issues?
mean(subset(le, le.group==1)$le.response) ## control group avg
sum(subset(le, le.group==1)$le.response==0)/nrow(le)
sum(subset(le, le.group==1)$le.response==4)/nrow(le)
sum(subset(le, le.group==2)$le.response==4)/nrow(le)
sum(subset(le, le.group==3)$le.response==4)/nrow(le)

#Model results 1 appendix B
## Basic list analysis
basic <- lm(le.response ~ le.group, dat=le)
basic.any <- lm(le.response ~ any, dat=le)

#Main text: figure 1
## Save as figures/list-experiment/Basic-regplot.pdf
dwplot(basic, 
       vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>%
  relabel_predictors(c(le.group2="Offered Bribe",
                       le.group3="Voted As Bribed")) +
  xlab("Proportion Reporting More Behaviors than the Control Condition")


## Control for demographics (randomization check, variance control)
##

###Model 2 appendix B model
summary(balance <- lm(le.response ~ le.group + region + gender + age + caste + religion + edu2 + rural.urban + polpart + wealth, 
           dat=le))

summary(lm(le.response ~ any + region + gender + age + caste + religion + readletter + writeletter + edu + house.material + rural.urban + cook, 
           dat=le))

#cited stat footnote 14 main text
## Correlation between wealth and education
cor(le$wealth, le$edu2)

## Subgroups
basic.women <- lm(le.response ~ le.group, dat=subset(le, gender=="g-female"))
basic.men <- lm(le.response ~ le.group, dat=subset(le, gender=="g-male"))
interact.gender <- lm(le.response ~ le.group*gender, dat=le)

basic.secondary <- lm(le.response ~le.group, dat=subset(le, edu2==TRUE))
basic.ltsecondary <- lm(le.response ~le.group, dat=subset(le, edu2==FALSE))
interact.edu <- lm(le.response ~ le.group*edu2, dat=le)

basic.urban <- lm(le.response ~le.group, dat=subset(le, rural.urban=="urban"))
basic.rural <- lm(le.response ~le.group, dat=subset(le, rural.urban=="rural"))
interact.urban <- lm(le.response ~ le.group*rural.urban, dat=le)

basic.polparthigh <- lm(le.response ~ le.group, dat=subset(le, polpart.high==TRUE))
basic.polpartlow <- lm(le.response ~ le.group, dat=subset(le, polpart.high==FALSE))
interact.polpart <- lm(le.response ~ le.group*polpart, dat=le)

basic.wealthhigh <- lm(le.response ~ le.group, dat=subset(le, wealth.high==TRUE))
basic.wealthlow <- lm(le.response ~ le.group, dat=subset(le, wealth.high==FALSE))
interact.wealth <- lm(le.response ~ le.group*wealth, dat=le)

interact.selfwealth <- lm(le.response ~ le.group*self.relwealth, dat=le)

#Model 3 appendix
## Model for statsig in appendix
interact.all <- lm(le.response ~ le.group*gender + le.group*edu2 + le.group*rural.urban + le.group*polpart + le.group*wealth, data=le)
stargazer(basic, balance, interact.all)

summary(interact.all)

le$ngo <- FALSE
le$ngo[grep("non-governmental", le$who.conduct)] <- TRUE
interact.ngo <- lm(le.response ~ le.group*ngo, dat=le)

le$io <- FALSE
le$io[grep("international", le$who.conduct)] <- TRUE
interact.io <- lm(le.response ~ le.group*io, dat=le)

subgroup.models <- bind_rows(map(list(basic.men, basic.women,
                                      basic.ltsecondary, basic.secondary,
                                      basic.rural, basic.urban, basic.wealthlow,
                                      basic.wealthhigh, basic.polpartlow,
                                      basic.polparthigh), 
                             tidy), .id="model") 
subgroup.models %<>% mutate(model=recode(model,
                                 `1` = "Men",
                                 `2` = "Women",
                                 `3` = "<Secondary Ed.",
                                 `4` = ">=Secondary Ed.",
                                 `5` = "Rural",
                                 `6` = "Urban",
                                 `7` = "Low Wealth",
                                 `8` = "High Wealth",
                                 `9` = "Low Pol. Part.",
                                 `10` = "High Pol. Part."))

#Figure 3 main text
#pdf("newfig.pdf", width=8, height=4)
dwplot(subgroup.models, dot_args = list(aes(shape = model)))  %>% 
  relabel_predictors(c(le.group2="Offered Bribe",
                       le.group3="Voted As Bribed")) +
  xlab("Proportion Reporting More Behaviors than the Control Condition") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  scale_shape_manual(values=1:10) +
  guides(
    shape = guide_legend("Model", reverse=TRUE), 
    colour = guide_legend("Model", reverse=TRUE))

### Are the differences in the subgroup ed coefficients statsig?
## Model 1 difference
subgroup.models[8,3] - subgroup.models[11,3]
## SE of model 1 difference
sqrt(subgroup.models[8,4]^2 + subgroup.models[11,4]^2)
## Point estimate minus 1.96*SE > 0
subgroup.models[8,3] - subgroup.models[11,3] - 1.96*sqrt(subgroup.models[8,4]^2 + subgroup.models[11,4]^2)

## Model 2 difference
subgroup.models[9,3] - subgroup.models[12,3]
## SE of model 1 difference
sqrt(subgroup.models[9,4]^2 + subgroup.models[12,4]^2)
## Point estimate minus 1.96*SE > 0
subgroup.models[9,3] - subgroup.models[12,3] - 1.96*sqrt(subgroup.models[9,4]^2 + subgroup.models[12,4]^2)


### Descriptives on bribes offered in constituency
bribes <- dat$group.1.clientism.measures.9.cli.mea.3.cli.mea.bribe
brb <- data.frame(loan.access=grepl("access-to-loans", bribes),
                  alcohol=grepl("alcohol", bribes),
                  building.materials=grepl("building-materials", bribes),
                  cash=grepl("cash", bribes),
                  consumer.goods=grepl("consumer-goods", bribes),
                  food=grepl("food", bribes),
                  jobs=grepl("jobs", bribes))
apply(brb, 2, mean)

#Figure 4 main text
## Save as figures/list-experiment/bribes-offered.pdf
barplot(apply(brb, 2, mean), names.arg=c("Loan", "Alcohol", "Build Mat.", "Cash", "Cons. Goods",
        "Food", "Jobs"))


### Direct clientelism questions and attitudes
direct <- dat[,c("group.1.clientism.measures.9.cli.mea.1a",
                 "group.1.clientism.measures.9.cli.mea.1b",
                 "group.1.clientism.measures.9.cli.mea.1c",
                 "group.1.clientism.measures.9.cli.mea.sel1.cli.mea.2",
                 "group.1.clientism.measures.9.cli.mea.sel1.cli.mea.3",
                 "group.1.clientism.measures.9.cli.mea.sel1.cli.mea.4",
                 "group.1.clientism.measures.9.cli.mea.4.1234",
                 "group.1.clientism.measures.9.cli.mea.4.2413",
                 "group.1.clientism.measures.9.cli.mea.4.4321",
                 "group.1.clientism.measures.9.cli.mea.4.3142",
                 "group.1.clientism.measures.9.pol.att.2")]
names(direct) <- c("stones.yes", "stones.no", "stones.notmoved", "accept.bribe",
                   "other.affect", "rejpct.to.stop", "statement.rand1", 
                   "statement.rand2", "statement.rand3", "statement.rand4", "attitude")
#Part C Appendix
## Stone experiement
## conservative
summary(direct[direct$stones.notmoved==0,]$stones.yes/10)
## lenient
summary(direct$stones.yes/10)
## midling
summary(direct$stones.yes/(direct$stones.yes + direct$stones.no))

## Direct would you accept bribe q
summary(direct$accept.bribe)
summary(direct$accept.bribe) / nrow(direct)
1-sum(direct$accept.bribe %in% c("very-likely", "slightly-likely")) / sum(direct$accept.bribe %in% c("very-unlikely", "slightly-unlikely"))

## Would others refusing affect your accepting a bribe?
summary(direct$other.affect) / nrow(direct)

## What percentage of people would need to turn down bribes to stop it?
summary(direct$rejpct.to.stop)
summary(direct$rejpct.to.stop) / nrow(direct)

## Which of these three statements is closest to your own opinion? 
## Statement 1: Because it helps to meet the needs of citizens in the community, it is acceptable for politicians to offer things to voters in exchange for their vote.
## Statement 2: Because they are not yet public servants and are operating as free citizens, it is acceptable for politicians to offer things to voters in exchange for their vote.
## Statement 3: Because spending resources on voters today means fewer resources to provide things to the community in the future, it is not acceptable for politicians to offer things to voters in exchange for their vote. 
## Statement 4: Because it goes against the principles of democracy, it is not acceptable for politicians to offer things to voters in exchange for their vote. 

## Q allowed multi-response, which makes this hard to use


## A candidate for national parliament is offering people who pledge their votes [a small amount of money/assistance obtaining small loans]. This behavior is:
summary(direct$attitude)
summary(direct$attitude) / nrow(direct)
